PATTERN  SEARCH  METHODS  FOR  LINEARLY  CONSTRAINED 
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Abstract.  This  paper  deals  with  generalized  pattern  search  (GPS)  algorithms  for  linearly 
constrained  optimization.  At  each  iteration,  the  GPS  algorithm  generates  a  set  of  directions  that 
conforms  to  the  geometry  of  any  nearby  linear  constraints,  and  this  set  is  used  to  define  the  POLL  set 
for  that  iteration.  The  contribution  of  this  paper  is  to  provide  a  detailed  algorithm  for  constructing 
the  set  of  directions  at  a  current  iterate  whether  or  not  the  constraints  are  degenerate.  The  main 
difficulty  in  the  degenerate  case  is  in  classifying  constraints  as  redundant  and  nonredundant.  We  give 
a  short  survey  of  the  main  definitions  and  methods  concerning  redundancy  and  propose  an  approach, 
which  may  be  useful  for  other  active  set  algorithms,  to  identify  the  nonredundant  constraints. 
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1.  Introduction.  This  paper  continues  the  development  of  the  generalized  pat¬ 
tern  search  (GPS)  algorithms  [T][9]  for  the  linearly  constrained  optimization  problems 


min  f(x) , 

x£il 


(1.1) 


where  /  :  3?"  — >  3?  U  {oo}  may  be  discontinuous,  and  the  feasible  region  is  given  by 


o  =  {ser 


aj x  <  bj 


i  £  1}  =  {x  £  5ft"  :  ATx  <  b}, 


(1.2) 


where  a*  £  5ft”,  bi  £  5ft,  and  A  £  Qnx  lJl  is  a  rational  matrix.  In  mi,  the  feasible 
region  Q  is  defined  as  tt  =  {x  £  5ft”  :  l  <  Ax  <  it}  where  A  £  Qmxn  is  a  rational 
matrix,  l,u  £  {5ft  U  {±oo}}m,  and  l  <  u.  As  is  evident,  (1.2)  reduces  to  the  definition 
in  mm,  where  the  rth  row  ar  of  the  matrix  A  is  equal  to  some  ith  row  af  of  the 
matrix  AT  with  a  coefficient  of  +1  or  —1,  and  bi  =  ur  or  bi  =  —lr. 

We  target  the  case  when  the  function  f(x)  can  be  an  expensive  black  box,  or 
provides  few  correct  digits  and  may  fail  to  return  a  value  even  for  feasible  points 
x  £  Cl.  In  this  situation,  the  accurate  approximation  of  derivatives  is  not  likely  to  be 
practical.  The  GPS  algorithms  rely  on  decrease  in  f(x);  an  iterate  x^+i  £  with 
f{xk+ i)  <  f(xk)  is  considered  successful. 

Lewis  and  Torczon  [5]  introduced  and  analyzed  the  generalized  pattern  search 
for  linearly  constrained  minimization  problems.  They  proved  that  if  the  objective 
function  is  continuously  differentiable  and  if  the  set  of  directions  that  defines  a  local 
search  is  chosen  properly  with  respect  to  the  geometry  of  the  boundary  of  the  feasible 
region,  then  the  GPS  converges  globally  to  a  Karush-Kuhn- Tucker  point.  Audet  and 
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Dennis  [T]  simplified  the  analysis  in  [9]  and  provided  new  convergence  results  for  a 
locally  Lipschitz  objective  function  by  applying  the  Clarke  generalized  directional 
derivative  [7]  to  the  pattern  search  methods. 

Generalized  pattern  search  algorithms  generate  a  sequence  of  iterates  {xk}  in  SR" 
with  non-increasing  objective  function  values.  In  linearly  constrained  optimization,  a 
set  of  directions  that  defines  the  so-called  poll  step  must  conform  to  the  geometry  of 
the  boundary  of  the  feasible  region.  The  key  idea,  which  was  first  suggested  by  May 
in  m  and  applied  to  the  GPS  in  [5],  is  to  use  as  search  directions  the  generators 
of  cones  that  are  polar  to  cones  generated  by  the  normals  of  faces  near  the  current 
iterate. 

Lewis  and  Torczon  [9j  showed  that  in  general  it  is  possible  to  find  the  set  of  gener¬ 
ators  by  vertex  enumeration  techniques  [2j,  but,  as  they  mentioned,  an  application  of 
this  approach  to  constructing  a  set  of  generators  can  be  expensive  in  the  degenerate 
case  when  the  active  constraints  are  linearly  dependent.  Lewis  and  Torczon  [S]  gave 
an  algorithm,  based  on  the  LU  factorization,  for  constructing  the  set  of  generators  in 
the  nondegenerate  case,  and  left  the  degenerate  case  for  future  work. 

Price  and  Coope  [12]  gave  as  an  aside  a  result  that  can  be  used  for  constructing  a 
set  of  generators  in  the  degenerate  case.  It  follows  from  their  result  that,  in  order  to 
construct  a  set  of  generators,  it  is  sufficient  to  consider  maximal  linearly  independent 
subsets  of  the  active  constraints.  However,  this  approach  implies  enumeration  of  all 
possible  linearly  independent  subsets  of  maximal  rank  and  does  not  take  into  account 
properties  of  the  problem  that  can  help  to  reduce  this  enumeration.  Price  and  Coope 
|12lj  outlined  an  algorithm  for  constructing  frames,  but  it  was  not  their  point  to  work 
out  details  of  the  numerical  implementation  in  the  degenerate  case. 

The  purpose  of  this  paper  is  to  give  detailed  consideration  to  the  GPS  in  the 
degenerate  case.  Our  purpose  here  is  complementary  to  [T]  and  [9].  Our  main  result 
is  a  detailed  algorithm  for  constructing  the  set  of  generators  at  a  current  GPS  iterate 
in  both  the  degenerate  and  nondegenerate  cases.  To  construct  the  set  of  generators  in 
the  degenerate  case,  we  identify  the  redundant  and  nonredundant  active  constraints 
and  then  use  either  QR  or  LU  decomposition. 

Classification  of  constraints  as  redundant  or  nonredundant  is  one  of  the  main 
issues  here,  because  it  is  sufficient  to  construct  the  set  of  generators  only  for  nonre¬ 
dundant  constraints.  Several  methods  to  classify  constraints  exist.  For  example,  there 
are  deterministic  algorithms  probabilistic  hit-and-run  methods  [3] ,  a  probabilis¬ 
tic  method  based  on  an  equivalence  between  the  constraint  classification  problem  and 
the  problem  of  finding  a  feasible  solution  to  a  set  covering  problem  [5] .  A  survey  and 
comparison  of  strategies  for  classifying  constraints  are  given  in  M-  Any  of  these 
approaches  can  be  applied  in  the  framework  of  the  GPS  to  identify  redundant  and 
nonredundant  constraints.  However,  in  the  paper,  we  propose  a  new  projection  ap¬ 
proach  to  identify  nonredundant  constraints,  that  is  more  suitable  for  GPS  methods. 

The  projection  method  is  similar  to  the  hit-and-run  algorithm  [2]  in  which  nonre¬ 
dundant  constraints  are  searched  for  along  random  direction  vectors  from  each  point 
in  a  sequence  of  random  interior  points.  In  contrast  to  the  hit-and-run  algorithm,  the 
projection  method  searches  for  a  nonredundant  constraint  in  a  deterministic  direction. 

The  major  advantage  of  the  projection  method  for  our  application  is  that  the 
number  of  direction  vectors  (in  the  terminology  of  the  hit-and-run  algorithm)  is  equal 
to  the  number  of  constraints  that  have  to  be  identified.  For  us  this  is  generally  a 
small  number.  In  the  hit-and-run  algorithm,  this  number  is  determined  by  a  stop 
criterion  and  can  be  large  if  many  of  the  random  generated  directions  do  not  detect  a 
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nonredundant  constraint.  Moreover,  the  formulas  used  in  the  projection  method  are 
simpler  than  those  used  for  computing  the  intersection  points  of  a  direction  vector 
with  the  hyperplanes  in  the  hit-and-run  algorithm.  Let  us  note  that  the  goal  of  the  hit- 
and-run  algorithm  is  detecting  all  nonredundant  constraints  in  a  full  system  of  linear 
inequalities.  We  use  the  projection  method  to  detect  the  nonredundant  constraints 
among  only  active  constraints  in  the  case  when  they  are  linearly  dependent. 

As  our  numerical  tests  show,  the  projection  method  cheaply  detects  all,  or  almost 
all,  nonredundant  constraints.  To  classify  constraints  that  have  not  been  detected  by 
the  projection  method,  we  use  another  approach  outlined  in  [6].  As  a  result,  we  guar¬ 
antee  that  every  active  constraint  is  detected  as  either  redundant  or  nonredundant. 
In  the  worst  case,  some  vertex  enumeration  techniques  [2]  mentioned  in  [5]  might  be 
necessary,  but  our  procedure  for  classification  of  the  constraints  seems  to  eliminate 
this  expense  for  many  cases. 

The  organization  of  the  paper  is  as  follows:  in  the  next  section,  we  give  a  brief 
description  of  GPS  as  well  as  the  convergence  result  for  linearly  constrained  mini¬ 
mization  following  papers  by  Audet  and  Dennis  [T] ,  and  by  Lewis  and  Torczon  [9] . 

Section  [3]  is  devoted  to  the  topic  of  redundancy.  In  the  first  part  of  the  section, 
we  introduce  a  definition  of  the  e-active  constraints  and  discuss  some  scaling  issues. 
The  second  part  of  Section  [3]  contains  essential  definitions  and  results  concerning  re¬ 
dundancy  0 1  m  ng  required  for  our  analysis.  Then  we  propose  our  projection 
method  to  determine  nonredundant  constraints,  and  we  briefly  describe  a  more  ex¬ 
pensive  follow  up  approach  to  be  applied  if  some  constraints  are  not  identified  by  the 
projection  method. 

In  Section  |4j  we  give  an  algorithm  for  the  set  of  generators,  and  we  discuss 
implementation  details,  including  rationality.  Section[5]is  devoted  to  some  concluding 
remarks. 

2.  Generalized  pattern  search  algorithms.  In  this  section,  we  give  a  brief 
description  with  the  convergence  result  for  the  GPS  methods  for  linearly  constrained 
minimization.  We  follow  papers  by  Audet  and  Dennis  [I],  and  by  Lewis  and  Tor¬ 
czon  [5],  and  we  refer  the  reader  there  for  details  of  managing  the  mesh  size  A  j,. 
Throughout,  we  will  always  use  the  1 2  norm. 

The  GPS  algorithms  can  be  applied  either  to  the  objective  function  /  or  to  the 
barrier  function  /a  =  /  +  ifn  :  5ft  — >  3?  U  {+00},  where  ipci  is  the  indicator  function 
for  Cl,  which  is  zero  on  Cl  and  00  elsewhere.  The  value  of  /a  is  +00  on  all  points  that 
are  either  infeasible  or  at  which  /  is  declared  to  be  +00.  This  barrier  approach  is 
probably  as  old  as  direct  search  methods  themselves. 

A  GPS  algorithm  for  linearly  constrained  optimization  generates  a  sequence  of 
iterates  {xk}  in  Cl.  The  current  iterate  Xk  £  3?"  is  chosen  from  a  finite  number  of 
points  on  a  mesh ,  which  is  a  discrete  subset  of  3?".  At  iteration  k,  the  mesh  is  centered 
around  the  current  mesh  point  (current  iterate)  Xk  and  its  fineness  is  parameterized 
through  the  mesh  size  parameter  Ak  >  0  as  follows 

Mk  =  {xk  +  A kDz  :  z  £  Z^},  (2.1) 

where  z\ is  the  set  of  nonnegative  integers,  and  D  is  a  set  of  positive  spanning 
directions  in  3ffn.  At  each  iteration,  some  positive  spanning  matrix  D ^  composed  of 
columns  of  D  is  used  to  construct  the  poll  set  P&,  i.e., 

Pfc  —  {Afc  T  Afcd  .  d  £  Dk }  • 
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(2.2) 


Xk+Akd 

d(zDk 


If  Xk  £  is  not  near  the  boundary,  then  D is  a  positive  spanning  set  for  3?n 

If  Xk  £  is  near  the  boundary,  the  matrix  D k  is  constructed  so  its  columns  dj 
also  span  the  cone  of  feasible  directions  at  Xk  and  conform  to  the  geometry  of  the 
boundary  of  f l.  Hence,  the  set  D  must  be  rich  enough  to  contain  generators  for  the 
tangent  cone  Tq(x)  —  cl {/i(w  —  x)  :  /i  >  0,  u>  £  H}  for  every  x  £  fi,  as  in: 

Definition  2.1.  A  rule  for  selecting  the  positive  spanning  sets  Dk  C  D  conforms 
to  f l  for  some  e  >  0,  if  at  each  iteration  k  and  for  each  y  in  the  boundary  of  fi  for 
which  || y  —  Xk\\  <  £,  Tn(y)  is  generated  by  a  nonnegative  linear  combination  of  the 
columns  of  a  subset  D^  of  Dk- 

Each  GPS  iteration  is  divided  into  two  phases:  an  optional  search  and  a  local  poll. 
In  the  SEARCH  step,  the  barrier  objective  function  /q  is  evaluated  at  a  finite  number 
of  points  on  a  mesh  to  try  to  find  one  that  yields  a  lower  objective  function  value  than 
the  incumbent.  This  step  is  extremely  flexible,  and  it  can  be  used  to  incorporate  the 
user’s  domain  knowledge. 

When  the  incumbent  is  replaced,  i.e. ,  when  fn(xk+i)  <  fn(%k),  or  equivalently 
when  f(xk+ i)  <  f(xk),  then  Xk+i  is  said  to  be  an  improved  mesh  point.  When  the 
search  step  fails  in  providing  an  improved  mesh  point,  the  POLL  step  is  invoked.  This 
second  step  consists  of  evaluating  the  barrier  objective  function  at  the  neighboring 
mesh  points  to  see  if  a  lower  function  value  can  be  found  there.  When  the  POLL  step 
fails  in  providing  an  improved  mesh  point,  then  the  current  incumbent  solution  is  said 
to  be  a  mesh  local  optimizer. 

We  remind  the  reader  that  the  normal  cone  Nq(x)  to  fl  at  x  is  the  nonnegative 
span  of  all  the  outwardly  pointing  constraint  normals  at  x  and  can  be  written  as  the 
polar  of  the  tangent  cone:  Nq(x)  =  {u  £  3?ra  :  Vw  £  Tq(x):  vtlu  <  0}. 

Assumptions.  We  make  the  following  standard  assumptions  [T] : 

Al:  A  function  fn  and  Xq  £  3?n  (with  fu(xo)  <  oo)  are  available. 

A2:  The  constraint  matrix  A  is  rational. 

A3:  All  iterates  {a;*,}  produced  by  the  Generalized  Pattern  Search  (GPS)  algorithm 
lie  in  a  compact  set. 

Audet  and  Dennis  [T]  proved  the  following  convergence  result  for  the  GPS  in  the 
linearly  constrained  case  using  only  these  assumptions. 

Theorem  2.2  (convergence  to  a  Karush-Kuhn-Tucker  point).  (Tj  Under  as¬ 
sumptions  A1-A3,  if  f  is  strictly  differentiable  at  a  limit  point  x  of  a  subsequence  of 
{xk}  at  which  is  decreased  and  for  which  the  corresponding  subsequence  of  {A^} 
goes  to  0,  and  if  the  rule  for  selecting  the  positive  spanning  sets  Dk  C  D  conforms  to 
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Initialization: 

Let  Xo  be  such  that  /n(^o)  is  finite.  Let  D  be  a  positive  spanning  set,  and 
let  M0  be  the  mesh  on  3?”  defined  by  A0  >  0,  and  D0.  Set  the  iteration 
counter  k  to  0. 

Search  and  poll  step: 

Perform  the  SEARCH  and  possibly  the  poll  steps  (or  only  part  of  them) 
until  an  improved  mesh  point  xk+\  with  the  lowest  /n  value  so  far  is  found 
on  the  mesh  M k  defined  by  equation  (2.1 1. 

—  Optional  SEARCH:  Evaluate  fa  on  a  finite  subset  of  trial  points  on  the 
mesh  Mk  defined  by  equation  (2.1)  (the  strategy  that  gives  the  set  of 
points  is  usually  provided  by  the  user;  it  must  be  finite  and  the  set 
can  be  empty). 

—  Local  POLL:  Evaluate  fa  on  the  poll  set  defined  in  equation  (2.2). 
Parameter  update: 

If  the  SEARCH  or  the  poll  step  produced  an  improved  mesh  point,  i.e.,  a 
feasible  iterate  x^+i  £  M k  n  for  which  fa(xk+i)  <  fa(xk),  then  update 


Afc+i  ■">  A^. 

Otherwise,  fa(xk )  <  fa{xk  +  Afcd)  for  all  d  £  Dk  and  so  xk  is  a  mesh  local 
optimizer.  Set  xk+i  =  xk,  update  Afc+i  <  A^. 

Increase  k  <—  k  +  1  and  go  back  to  the  SEARCH  and  poll  step. 


Fig.  2.1.  A  simple  GPS  algorithm 


Q  for  some  e  >  0,  then  V f(x)Tu>  >  0  for  all  to  £  Tq(x)  and  so  — V/(x)  £  Nq(x). 
Thus,  x  is  a  Karush-Kuhn-  Tucker  point. 

The  purpose  of  this  paper  is  to  provide  an  algorithm  for  constructing  sets  Dk 
that  conform  to  the  boundary  of  fi.  If  the  active  constraints  are  linearly  dependent, 
we  apply  strategies  for  the  identification  of  redundant  and  nonredundant  constraints, 
which  are  described  in  the  next  section,  and  then  construct  sets  Dk  taking  into  account 
only  nonredundant  constraints.  We  now  pause  to  outline  the  main  results  concerning 
redundancy  from  mathematical  programming,  and  then  in  Section  [4]  we  continue 
consideration  of  the  GPS  and  strategies  for  constructing  the  sets  Dk. 


3.  Redundancy.  In  this  section,  we  give  essential  definitions  and  results  con¬ 
cerning  redundancy  ®  ©  0  HU  ng  that  are  required  for  our  analysis.  Then  we 
propose  our  approach,  the  projection  method,  to  determining  the  nonredundant  con¬ 
straints  and  briefly  describe  another  approach  that  is  applied  if  some  constraints  are 
not  identified  by  the  projection  method. 


We  consider  the  feasible  region  f l  defined  by  (1.2 1,  and  refer  to  the  inequality 
aj x  <  bj  as  the  j-th  constraint.  The  region  represented  by  all  but  the  j'th  constraint 
is  given  by 


flj  =  {x  £  3ffra  :  ajx  <  bi ,  i  £  I\j}, 

where  I\j  is  the  set  I  with  element  j  removed. 

3.1.  e— active  constraints.  In  this  subsection,  we  introduce  a  definition  of  the 
e-active  constraints  and  discuss  some  scaling  issues.  First  we  give  the  definitions  used 
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by  Lewis  and  Torczon  [9]  and  by  Price  and  Coope  [T2], 

Definition  3.1.  ( e.g .,  [12]).  Let  some  scalar  e  >  0  be  given  and  x k  £  ft.  The 
jth  constraint  is  e-active  at  Xk  if 

0  <  bj  —  ajxk  <  e.  (3.1) 

Definition  3.2.  (e.g.,  i)-  Let  some  scalar  e  >  0  be  given  and  Xk  £  LI.  The  jth 
constraint  is  e-active  at  Xk  if 


dist(£fc,  Hj)  <  e, 


(3.2) 


where  Hj  =  {x  £  Sftn  :  ajx  =  bj},  and  dist (xk,Hj)  =  min  \\y  —  Xk\\  is  the  distance 
from  Xk  to  the  hyperplane  Hj . 

As  is  easy  to  see,  the  jth  constraint  can  be  made  e-active  at  Xk  in  the  sense 
of  Definition  3.1  by  multiplying  the  inequality  bj  —  ajxk  >  0  by  a  sufficiently  small 
number.  On  the  other  hand,  this  multiplication  does  not  change  the  distance  between 
the  point  Xk  and  any  Hj  defined  in  Definition  3.2  In  the  paper,  we  prefer  to  use 


Definition  3T  since  it  is  easier  to  check  than  Definition  T2  However,  Definition  |3.1| 
is  proper,  if  we  assume  preliminary  scaling  of  the  constraints  so  that  the  following 
lemma  applies. 

Lemma  3.3.  Let  some  scalar  e  >  0  be  given,  Xk  £  LI,  and  1 1 a, y  1 1  =  1  for  all  j  £  I 
in  ( 1.2 1 .  Then,  for  any  j  £  I,  Definition  3.1  of  the  e-active  constraint  is  equivalent 


to  Definition  3.2  and  the  projection  Pj(xk)  of  the  point  Xk  onto  the  hyperplane  Hj  = 


{i  €  :  ajx  =  bj}  is  defined  by 


Pj(xk)  =  xk  +  aj(bj  -  a j  xk). 


(3.3) 


Proof.  For  any  j  £  I,  the  distance  from  Xk  to  the  hyperplane  Hj  =  {x  £  : 


ajx  =  bj }  is  given  by 


dist  (xk,Hj)  = 


I  bj  -  ajxk  | 


(3.4) 


Hence,  if  1 1 a,y 1 1  =  1  and  Xk  £  Ll,  (3.1 1  is  equivalent  to  (3.2 1. 

By  definition  of  the  projection  of  Xk  onto  Hj, 

\\Pj(xk)  -  xk\\  =  dist (xk,Hj). 

Since  Xk  £  Ll  and  1 1 a. y  1 1  =  1,  it  follows  from  (|3.4|)  that  dist(xfe,  Hj)  =  bj  —  ajxk  and 
Pj(xk)  =  xk  +  aj  dist(*fc,  Hj)  =  xk  +  aj(bj  -  ajxk). 


Hence,  (3.3)  holds.  □ 


To  satisfy  the  conditions  of  Lemma  3.3  we  introduce  the  matrix  A  that  is  an 


additional  scaled  copy  of  the  matrix  A  given  in  (1.2 1  and  a  scaled  vector  b  such  that 

b. i 


b,= 


at 


i  £  I. 


(3.5) 


Consequently,  ||aj||  =  1  for  all  i  £  I  and  Ll  =  {x  £  3?"  :  ATx  <  6}  =  {x  £  3?" 
AT x  <  b}  =  {s  £  3?"  :  ajx  <bi,  i  £  I}. 
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We  use  the  matrix  A  and  the  vector  b  to  define  the  set  of  indices  of  the  £-active 
constraints: 

I(xk,e)  =  {ie  I  :  0  <bi-  af  xk  <  e},  (3.6) 

as  well  as  to  apply  the  projection  method  for  detection  of  the  nonredundant  constraints 


(see  Section  3.3.1  for  more  details.)  We  refer  to  the  set  I(xk,£ )  as  the  working  index 
set  at  the  current  iterate  xk- 

The  paper  also  makes  use  of  the  regions  given  by 

f l(xk,e)  =  {a;  G  :  ajx  <bt,  i  €  I(xk,£)},  (3.7) 


flj(xk,e)  =  {x  G  3T  :  af  x  <bi:  i  G  I(xk,e)\j},  j  G  I{xk,e). 

Clearly,  f l  C  fl(xk,e)  C  flj(x :k,e). 

3.2.  Redundancy  in  mathematical  programming.  In  this  subsection,  we 
give  definitions  and  theorems  consistent  with  the  mathematical  programming  litera¬ 
ture  mm  mug. 

The  following  definitions  and  results  are  from  Dung. 

Definition  3.4  (Polyhedron).  A  subset  offtn  described  by  a  finite  set  of  linear 
constraints  P  =  {x  G  3Jn  :  Cx  <  d}  is  a  polyhedron. 


Obviously,  f l  given  by  (1.2 1  and  H(xk,e)  given  by  (|3.7|)  are  polyhedra. 
Definition  3.5.  The  points  x1, 


,  xk  G  3?"  are  affinely  independent  if  the  k  —  1 
directions  x2  —  x1, . . . ,  xk  —  x1  are  linearly  independent,  or  alternatively  the  k  vectors 
[x1, 1), . . . ,  (xk,  1)  G  3?"+1  are  linearly  independent. 

Definition  3.6.  The  dimension  of  P,  denoted  dim(P),  is  one  less  than  the 
maximum  number  of  affinely  independent  points  in  P . 

This  means  that  P  C  3f”  is  full- dimensional  if  and  only  if  dim(P)  =  n.  We  will 
assume  that  fl  is  full-dimensional.  If  it  were  not,  then  a  barrier  GPS  approach  would 
not  be  a  reasonable  way  to  handle  the  constraints  because  it  would  be  difficult  to  find 
SEARCH  or  POLL  points  in  Q.  Since  we  assume  is  full  dimensional,  this  implies  that 
its  supersets  n(xk,s)  and  Clj(xk,s)  are  full-dimensional. 

Definition  3.7  (Valid  inequality).  An  inequality  CjX  <  dj  is  a  valid  inequality 
for  P  C  if  CjX  <  dj  for  all  x  G  P. 

Definition  3.8  (Face  and  Facet),  (i)  P  defines  a  face  of  the  polyhedron  P  if 
F  =  {x  G  P  :  CjX  =  dj}  for  some  valid  inequality  CjX  <  dj  of  P.  F  0  is  said  to  be 
a  proper  face  of  P  if  F  P. 

(ii)  F  is  a  facet  of  P  if  F  is  a  face  of  P  and  dim(P)  =  dim(P)  —  1. 

Definition  3.9  (Dominance  of  inequalities).  If  CiX  <  di  and  CjX  <  dj  are  two 
valid  inequalities  for  P  C  3?" ,  C\X  <  di  dominates  CjX  <  d 
that  Ci  >  ucj  and  di  <  udj,  and  ( Ci,di )  ( ucj,udj ). 

Definition  3.10  (Redundant  inequality).  A  valid  inequality  CjX  <  dj  is  redun¬ 
dant  in  the  description  of  P  (in  other  words,  the  jth  constraint  is  redundant)  if  there 
exist  k  >  1  valid  inequalities  CiX  <  di  for  i  =  1, . . . ,  k  for  P,  and  weights  oti  >  0  for 
i  =  1, . . . ,  k  such  that  (X^=i  oiicfjx  <  Yli=i  dominates  CjX  <  dj. 

The  following  example  illustrates  that  we  can  not  replace  a.i  >  0  with  at  >  0  in 
Definition  13.101 

Example  1.  Let  the  following  inequalities  in  R2  be  valid  inequalities  for  some 
polyhedron  P: 


.jj.  ^  u,j  if  there  exists  u  >  0  such 


c\X  <  d\ 
c2x  <  d2 


X\<1 
-X\  <  1 


(3.8) 
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3.10 


If  ol\  were  equal  to  zero  in  Definition 
and 


we  would  obtain  a\C\  =  (0,0),  a±di  =  0, 


(0, 0)  =  aqci  >  UC2  =  (— 1,  0),  0  =  a\d\  <  ud,2  =  1,  u=  1. 


Hence,  by  Definition  3.9  the  first  inequality  in  (3.8 1  would  dominate  the  second 


one  and  we  would  make  the  wrong  conclusion,  by  Definition  (3.10),  that  the  second 
inequality  is  redundant  in  the  description  of  P.  □ 

Let  us  note  that  since  H  C  f2(a;fc,e),  if  the  jth  constraint  is  redundant  in  the 
description  of  fl(xk,e),  it  is  redundant  in  the  description  of  fl. 


Fig.  3.1.  An  illustration  of  e-active  and  redundant  constraints.  Constraints  1,  2,  and  3  are 
e-active  at  the  current  iterate  x  and  constraint  2  is  redundant. 


Definition  3.11  (Interior  point).  A  point  x  £  P  is  called  an  interior  point  of  P 
if  Cx  <  d. 

We  need  the  following  results  from  integer  programming  m  pp.  142-144]  and 
[Til  pp.  85-92]. 

Proposition  3.12.  m  Corollary  2.5]  A  polyhedron  is  full  dimensional  if  and 
only  if  it  has  an  interior  point. 

Theorem  3.13.  [15.,  Theorem  9.1]  If  P  is  a  full- dimensional  polyhedron,  it  has 
a  unique  minimal  description 


P  =  {x£$in 


CiX  <  di , 


i=  1,  •  •  • ,  m}, 


where  each  inequality  is  unique  to  within  a  positive  multiple. 

Corollary  3.14.  [15j  Proposition  9.2]  If  P  is  full- dimensional,  a  valid  inequality 
CjX  <  dj  is  necessary  in  the  description  of  P  if  and  only  if  it  defines  a  facet  of  P. 


fined  in  (3.7). 


Corollary  3.14  means  that  the  following  concepts  are  equivalent  for  H(xfc,e)  de- 


•  The  jth  inequality  ajx  <  bj  defines  a  facet  of  H(xfc,e). 

•  The  jth  inequality  ajx  <  bj  is  nonredundant  in  the  description  of  V,(xk,e). 

•  The  jth  inequality  ajx  <  bj  is  necessary  in  description  of  fl(xk,  e),  or  in  other 
words, 


&{xk,e)  C  nj(xk,e).  (3.9) 

Our  approach  to  the  identification  of  nonredundant  constraints  is  based  primarily 
on  the  following  proposition. 

Proposition  3.15.  Let  a  working  index  set  I(xk,e)  be  given.  An  inequality 
ajx  <  bj,  j  £  J(xfc,e),  is  nonredundant  in  the  description  ofH(xk,e)  if  and  only  if 


either  I(xk,e)  =  {j}  or  there  exists  x  G  R”  such  that  afx  =  bj  and  af  x  <  bi  for  all 
i  G  I(xk,e)\j. 

Proof.  Since  the  case  I(xk,  e)  =  {j}  is  trivial,  we  give  the  proof  for  the  case  when 
I(xk,e)\j  /  0. 

Necessity.  Since  the  inequality  af  x  <  bj  is  nonredundant,  then,  by  (3.9),  there 


exists  x*  G  Rn  such  that  afx*  <  bi  for  all  i  G  I(xk,s)\j,  and  afx*  >  bj.  By 


Proposition 


3.12 


there  exists  an  interior  point  x  G  £l(xk,s)  such  that  afx  <  bi  for 


all  i  G  I(xk,e).  Thus  on  the  line  between  x*  and  x  there  is  a  point  igF  satisfying 
afx  =  bj  and  afx  <  bi  for  all  i  G  I(xk,s)\j. 

Sufficiency.  To  obtain  a  contradiction,  we  suppose  that  there  exists  iGS"  such 
that  afx  =  bj  and  afx  <  bi  for  all  i  G  I(x.k,s)\j,  but  the  inequality  afx  <  bj 
is  redundant.  Then,  using  Definition  |3.10[  we  obtain  that  there  exist  k  >  1  valid 
inequalities  afx  <  bi  for  i  =  1, . . . ,  k  for  f l(xk,  e),  and  weights  ai  >  0  for  i  =  1, . . . ,  k 
such  that  iffi=iaiaf)x  —  dominates  afx  <  bj.  Thus,  by  Definition 

there  exists  u  >  0  such  that 


3.9 


uajX  <  ataf)x  <  <  ubj. 


i- 1 


i=l 


(3.10) 


Since  by  hypothesis,  afx  <  bi  for  all  i  j,  it  follows  from  (3.101  that 

k  k 


(^aiiaf) x  <  y^  ajbj. 


Hence,  by  ( 

3.10 

,  ajX 

Proposition 

3.15 

This  contradicts  the  supposition  that  ajX  =  bj.  □ 


then  there  exists  a  feasible  point  x  G  fl(xk,s)  such  that  only  this  constraint  holds 
with  equality  at  x. 

Our  approach  to  the  identification  of  the  redundant  constraints  is  based  primarily 
on  the  following  theorem  [5]. 

Theorem  3.16. 


The  jth  constraint  is  redundant  in  system  (1.2 1  if  and  only  if 

(3.11) 


maximize  aj  x, 


subject  to  x  G  flj. 


has  an  optimal  solution  x*  such  that  afx*  <  bj. 

3.3.  Approaches  to  identification  of  redundant  and  nonredundant  con¬ 
straints.  In  section  [4j  to  identify  the  redundant  constraints,  we  use  an  approach 
proposed  in  |6j,  based  on  Theorem  3.16  and  briefly  outlined  in  Section  3.3.2  and  to 
determine  the  nonredundnat  constraints,  we  apply  a  method  presented  in  the  next 
subsection. 

3.3.1.  A  projection  method.  In  this  subsection,  we  propose  the  projection 
method  that  is  intended  to  identify  nonredundant  constraints.  The  main  idea  of  this 
method  is  to  construct,  if  it  is  possible,  a  point  x  such  that  afx  =  bj  and  afx<bi  for 
all  i  G  I(xk ,  e)\j.  If  such  a  point  x  exists,  then  by  Proposition  3.15  the  jth  constraint 
is  nonredundant. 


We  recall  that,  in  (3.5),  we  defined  a  scaled  copy  A  of  the  matrix  A  and  a  scaled 
vector  b.  We  denote  by  Pj(xk ),  the  projection  of  the  point  Xk  G  SR"  onto  the  hy¬ 


perplane  Hj  =  {i  G  S"  :  afx  =  bj}.  Assume  that  Xk  G  D.  Then  by  (3.3 1  and  by 
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IMI  =  i 


Pj(xk)  =  xk  +  aj(bj  -  ajxk).  (3-12) 

The  following  proposition  is  the  main  one  for  the  projection  method. 
Proposition  3.17.  Let  xk  €  LI  and  let  a  working  index  set  I(xk,e)  be  given.  An 
inequality  ajx  <  bj,  j  £  I{xk,e),  is  nonredundant  in  the  description  of  Q(xk,e)  if 

aJPj{xk)<bi  for  all  i£l(xk,e)\j,  (3.13) 


where  Pj{xk )  is  a  projection  of  xk  onto  Hj. 

Proof.  The  proof  follows  from  Proposition  |3.15|  □ 

An  application  of  the  projection  method  follows  from  Proposition  jum  Namely, 
we  classify  the  jth  constraint  as  nonredundant,  where  j  £  I(xk,s),  if  (3.13)  holds  for 
all  i  £  I(xk,e)\j  with  P, 


;  j(xk)  given  by  (3.12). 


3.3.2.  Approach  for  identifying  redundant  and  nonredundant  const¬ 
raints.  If  some  constraints  have  not  been  identified  by  the  projection  method,  we 
apply  another  approach  based  on  Theorem  |3.16| to  identify  redundant  and  nonredun¬ 
dant  constraints.  The  book  [Hj  contains  description  of  different  methods  in  the  context 
of  this  approach.  These  methods  use  some  very  special  propositions  concerning  slack 
variables  and  allow  the  authors  to  simplify  numerical  solution  of  the  linear  program¬ 


ming  (LP)  problem  (3.11)  or  to  substitute  solving  the  LP  problem  by  some  other 
procedure,  which  is  less  expensive  from  the  numerical  point  of  view.  We  have  not 
included  these  details  of  numerical  implementation  of  the  solution  of  the  LP  problem 
into  the  paper,  since  the  reader  can  find  them  in  j8]. 

4.  Construction  of  the  set  of  generators.  The  purpose  of  this  section  is  to 
provide  a  detailed  algorithm  for  constructing  the  set  Dk  introduced  in  Section  [2] 

Let  some  scalar  e  >  0  be  given,  and  let  af  be  the  zth  row  of  the  matrix  AT  in 


(3.5 ).  At  the  current  iterate  xk,  we  construct  the  working  index  set  I(xk,  e)  as  follows: 


0  <  bi  —  aj xk  <  e 


i  £  I(xk,e). 


The  last  inequality  means  that  every  constraint  that  is  active  at  xk  or  at  some  point 
near  xk  appears  in  I(xk,e).  In  [I],  the  authors  suggest  not  setting  e  so  small  that  Ak 
is  made  small  by  approaching  the  boundary  too  closely  before  including  conforming 
directions  that  allow  the  iterates  to  move  along  the  boundary  of  f 2. 

Without  loss  of  generality,  we  assume  that  I(xk,  s)  =  {1, . . . ,  m},  for  m  >  2.  This 
avoids  more  cumbersome  notation,  like  I(xk,s)  =  {i i(xk,e), . . .  ,im(xk,£)}. 

We  denote  by  Bk,  the  matrix  whose  columns  are  the  columns  of  A  corresponding 
to  the  indices  I(xk,e)  =  {1, . . . ,  m};  i.e., 

Bk  =  [a1...am\.  (4.1) 


4.1.  Classification  of  degeneracy  at  the  current  iterate.  Let  the  matrix 
Bk  be  defined  by  (4.1 1.  At  the  current  iterate  xk,  one  of  the  following  cases  holds: 

•  a  nondegenerate  case  when  the  matrix  Bk  has  full  rank; 

•  a  degenerate  case  when  Bk  does  not  have  full  rank  and  all  nonredundant  constraints 
are  linearly  independent; 

•  a  degenerate  nonredundant  case  when  Bk  does  not  have  full  rank  and  all  nonredun¬ 
dant  constraints  are  linearly  dependent.  This  case  can  be  illustrated  by  the  following 
example. 
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Example  2.  (Charles  Audet  sent  this  example  to  us.)  Suppose  that  the  feasible 


region  Q  (1.2 1  is  defined  by  the  following  system  shown  in  Figure  4.1 


XI  —  2X2  —  2x3 

< 

0 

—2xi  +  x2  -  2x3 

< 

0 

—2xi  -  2x2  +  x3 

< 

0 

Xi 

> 

0 

x2 

> 

0 

X3 

> 

0 

If  Xk  G  3?3  is  near  the  origin,  all  six  constraints  are  active,  linearly  dependent, 
and  nonredundant.  The  matrix  B k  is  given  as 

/  1  -2  -2  -1  0  0  \ 

Bk  =  -2  1  -2  0  -1  0  . 

\  -2  -2  1  0  0  -1  / 


Fig.  4.1.  Example  2.  An  illustration  of  the  degenerate  nonredundant  case. 


4.2.  Set  of  generators.  Following  [SJ,  we  define  the  cone  K{x k,£)  as  the  cone 
generated  by  the  normals  to  the  e-active  constraints,  and  K°(xk,£ )  as  its  polar: 


K°(xk,e )  =  {w  G  $tn  :  a[ w  <0  V*  €  I(xk,£)}. 


(4.3) 


M- 


The  cone  K°(xk,e)  defined  by  (4.3)  can  be  rewritten  as  a  finitely  generated  cone 


K°(xk)e)  =  {w  :  w  =  ^2XjVj,  \j>0,  j  = 
i= i 


(4.4) 


where  the  vectors  v\, .  ..,vr  are  a  set  of  generators  for  the  cone  K°(xk,£)  defined  as 
follows. 

Definition  4.1  (Set  of  generators).  A  set  V  =  (ui, . . . ,  ur}  is  called  a  set  of 
generators  of  the  cone  K0(xk,£)  defined  by  (4.3 1  if 

(i)  any  vector  v  G  K°(xk,£)  can  be  represented  as  a  nonnegative  linear  combination 
of  vectors  v,;  from  V,  i.e.  (|4.4|)  holds; 

(ii)  no  subset  of  {v i, . . .  ,vr}  satisfies  (4.4 1 . 
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The  key  idea,  which  was  first  suggested  by  May  in  [ICt  and  applied  to  the  GPS 
in  [5],  is  to  include  in  Dk  the  generators  of  the  cone  K°(x fc,e).  Hence,  the  problem 
of  construction  of  the  set  Dk  reduces  to  the  problem  of  constructing  the  generators 
{i>i, . . . ,  vr}  of  the  cone  K°(xk,e)  and  then  completing  them  to  a  positive  spanning 
set  for  5ft". 

The  following  proposition  means  that  it  is  sufficient  to  construct  the  set  of  gen¬ 
erators  only  for  nonredundant  constraints. 

Proposition  4.2.  Let  I{xk,e)  be  a  set  of  indices  of  constraints  that  are  e-active 
at  the  point  xk.  Let  Iisr(xk,s)  C  I{xk,e)  be  the  subset  of  indices  of  the  nonredundant 
constraints  that  define  Ll(xk,e).  Let  the  cone  K0(xk,s)  be  defined  by  (4.3 1  and  let  the 
cone  K^{xk,e)  be  given  by 

K%(xk,e)  =  {ro  €  S"  :  af  w  <0  V*  e  IN(xk,e)}. 

If  {v i , . . . ,  vp}  is  a  set  of  generators  for  K{fr{xk,e),  then  the  set  of  vectors  {vi, . . . ,  vp} 
is  a  set  of  generators  for  K°(xk,e). 

Proof.  The  proof  of  this  proposition  follows  from  Corollary  |3.14|  □ 

As  is  mentioned  in  |9j,  pattern  search  methods  require  their  iterates  to  lie  on  a 
rational  lattice.  To  arrange  this,  Lewis  and  Torczon  [9]  placed  an  additional  require¬ 
ment  that  the  matrix  of  constraints  AT  in  is  rational.  Under  this  requirement, 
Lewis  and  Torczon  [9]  showed,  in  the  following  theorem,  that  it  is  always  possible  to 
find  rational  generators  for  the  cones  K°(xk,e),  which,  with  the  rational  mesh  size 
parameter  A*,,  makes  the  GPS  points  lie  on  a  rational  lattice. 

Theorem  4.3.  Suppose  K  is  a  cone  with  rational  generators  V.  Then  there 
exists  a  set  of  rational  generators  for  K° . 

Moreover,  for  the  case  of  linearly  independent  active  constraints,  Lewis  and 
Torczon  [5]  proposed  constructing  the  set  of  generators  for  all  the  cones  K(xk,e), 
0  <  e  <  <5,  as  follows: 

Theorem  4.4.  Suppose  that  for  some  5,  K(x,5)  has  a  linearly  independent  set 
of  rational  generators  V.  Let  N  be  a  rational  positive  basis  for  the  nullspace  ofVT. 

Then ,  for  any  £,  0  <  e  <  S,  a  set  of  rational  generators  for  K°{x,  e)  can  be  found 
among  the  columns  of  N,  V(VTV)~1 ,  and  —V(yTV)~1. 

As  is  shown  in  [5],  the  matrix  N  can  be  constructed  by  taking  columns  of  the 
matrices  ±(J—  V(VTV)~1VT). 


Let  us  recall  that  we  use  the  scaled  matrix  A  defined  in  (3.5 1  to  determine  £-active, 
redundant,  and  nonredundant  constraints.  Then  we  use  the  result  stated  in  Theorem 
|4.4|  together  with  rational  columns  of  A ,  which  correspond  to  the  nonredundant  and 
£-active  constraints,  to  obtain  a  set  of  rational  generators. 

A  set  of  generators,  which  may  be  irrational  in  exact  arithmetic,  can  also  be  found 
by  using  the  QR  factorization  of  the  matrix  V.  The  following  corollary  shows  how  to 
use  the  QR  factorization  of  V  to  construct  the  generators  for  all  the  cones  K°(xk,e), 
0  <  £  <  5.  We  recall  that  the  full  QR  factorization  of  V  can  be  represented  as 


V  —  [Q i  Q2] 


Ri 

0 


R-I 

0 


where  R\  is  upper  triangular  and  rank  R\  =  rank  V,  and  the  columns  of  Q±  form  an 
orthonormal  basis  for  the  space  spanned  by  the  columns  of  V,  while  the  columns  of 
Q 2  constitute  an  orthonormal  basis  for  nullspace  of  VT . 

Corollary  4.5.  Suppose  that  for  some  S,  K(x,5)  has  a  linearly  independent  set 
of  rational  generators  V.  Then,  for  any  e,  0  <  £  <  S,  a  set  of  generators  for  K°(x,e) 
can  be  found  among  the  columns  of  Q2,  Q\R\{Rf  Ri)^1 ,  and  —  QiRi(R.J . 
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Proof.  By  substituting  V  =  QR  and  using  the  properties  of  the  matrices  in  the 
QR  factorization,  we  obtain 


yT  dW 


VfV1  V)-1  =  QR{{QRy  ( QR ))'1  =  QR^R1  Q 1  QR)-1  =  QR{R 1  R) 


By  applying  Theorem  4.4  and  by  taking  into  account  that  columns  of  Q 2  span  the 
nullspace  of  VT ,  we  obtain  the  statement  of  the  corollary.  □ 

From  the  theoretical  point  of  view,  a  set  of  generators  obtained  by  using  Corollary 
4.5  may  be  irrational  since  an  implementation  of  the  QR  decomposition  involves 
calculation  of  square  roots.  Of  course,  we  use  floating  point  arithmetic,  and  so  these 
generators  are  rational,  but  they  probably  generate  a  different  cone  due  to  numerical 
errors.  Still  the  same  is  true  of  generators  found  by  using  the  LU  factorization. 

4.3.  An  algorithm  for  constructing  the  set  of  generators.  This  section 
presents  an  algorithm  for  constructing  a  set  of  generators  for  the  cone  K°(xk,£ )  at 
the  current  iterate  Xk  for  a  given  parameter  e. 

4.3.1.  Comments  on  the  algorithm.  The  algorithm  consists  of  two  main 
parts.  In  the  first  part,  we  determine  the  set  of  indices  of  the  nonredundant  e-active 
constraints  //v( Xk,£)  C  I(xk,s)  and  form  the  matrix  Bn  whose  columns  are  the 
columns  of  A  corresponding  to  the  indices  in  J/v  (£*,,£).  We  use  information  about 
the  set  I]y(xk,  s)  from  the  previous  iterations  of  the  GPS  algorithm.  Namely,  we  put 
into  the  set  J/v  (£&,£)  all  indices  that  correspond  to  the  e-active  constraints  at  the 
current  iterate  and  that  were  detected  as  indices  of  the  nonredundant  constraints  at 
the  previous  iterations  of  the  algorithm.  In  the  second  part  of  the  algorithm,  we 
construct  the  set  of  generators  Dk  required  by  the  GPS  and  by  Theorem  |2.2| 

First  of  all,  we  try  to  identify  the  nonredundant  active  constraints.  If  the  matrix 


Bk  defined  by  (4.1 1  has  full  rank,  then  all  e-active  constraints  are  nonredundant, 


IN(xk,e)  =  I(xk,e),  and  BN  =  Bk. 

If  the  matrix  Bk  does  not  have  full  rank  and  we  have  indices  that  have  not  been 
classified  at  the  previous  iterations  of  the  algorithm,  we  suggest  using  two  steps  in 
succession. 

The  first  strategy  is  intended  to  determine  nonredundant  constraints  cheaply  by 
applying  the  projection  method  described  in  section  nm  By  Proposition  |3.17j  if 
the  projection  Pj(xk)  of  the  current  iterate  Xk  onto  the  hyperplane  Hj  =  {x  €  JJn  : 
aj x  =  bj}  is  feasible,  and  only  the  jth  constraint  holds  with  equality  at  Pj(xk),  then 
the  jth  constraint  is  nonredundant,  and  we  can  put  index  j  into  the  set  /jv (£&,£). 
If  some  constraints  have  not  been  identified  by  the  projection  method,  we  can  either 
apply  the  projection  method  with  some  other  point  x  7^  Xk  or  apply  the  second 
strategy. 

The  second  strategy  is  intended  to  classify  redundant  and  nonredundant  con¬ 
straints  among  those  constraints  that  have  not  been  determined  as  nonredundant  by 
the  projection  method.  To  identify  the  constraint,  an  approach  outlined  in  j|6]  is  ap¬ 
plied.  The  basis  of  this  second  strategy  is  provided  by  Theorem  |3.16|  If  the  number 
of  constraints  we  have  to  identify  is  too  big,  we  can  skip  an  application  of  the  second 
strategy  to  the  algorithm  and  construct  a  set  of  generators  using  the  set  /jv(xfc,£) 
defined  by  application  of  the  first  strategy.  Then,  while  doing  the  poll  step,  if  we  find 
some  point  x  =  Xk  +  A  d,  where  d  is  some  column  of  Dk ,  such  that  ajx  >  bj  and 
afx  <  bi  for  all  i  £  I(xk,s)\j ,  we  can  conclude  that  fl(xk,£)  C  Clj(xk,£).  Hence,  by 


Corollary  3.14  the  jth  constraint  is  nonredundant  and  we  include  j  into  set  Ij\r(xk,s). 
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When  we  have  specified  all  redundant  and  nonredundant  constraints,  we  compose 
the  matrix  Bn  of  those  columns  of  the  rational  matrix  A  that  correspond  to  nonre¬ 
dundant  constraints.  To  determine  the  rank  of  Bn,  the  QR  factorization  can  be  used. 
If  Bn  has  full  rank,  we  can  apply  the  function  Set-QR  or  Set-LU,  which  uses  the 
QR  or  LU  decomposition  of  Bn,  to  construct  the  set  of  generators.  Both  functions 
are  constructed  on  the  basis  of  results  proved  in  Theorem  |4.4|  and  in  Corollary  |4.5| 
where  V  =  BN- 

If  the  matrix  Bn  does  not  have  full  rank,  we  apply  the  result  proved  in  m- 
Specifically,  to  construct  the  set  Dk,  it  is  sufficient  to  apply  the  function  Set_QR  or 
Set_LU  to  all  maximal  linearly  independent  subsets  of  the  columns  of  Bn-  We  can 
estimate  the  number  S  of  these  subsets.  If  Bn  has  rriN  columns  and  rank  r,  then 


S  = 


mj y! 

r!(mjv  —  r)! 


(4.5) 


In  this  worst  case,  some  vertex  enumeration  techniques  [2]  mentioned  in  [9]  might  be 
necessary,  but  our  procedure  for  classification  of  the  constraints  should  eliminate  this 
expense  for  many  cases. 


4.3.2.  Algorithm.  We  denote  the  set  of  indices  of  the  nonredundant  active 
constraints  by  lN{xk,£).  Thus,  for  j  £  I(xk,e), 

(i)  if  j  £  lN(xk,s),  the  inequality  ajx  <  bj  is  nonredundant;  and 

(ii)  if  j  £  I(xk,s)\lN{xk,£),  the  inequality  ajx  <  bj  is  redundant. 

We  use  In  Q  I  to  denote  the  set  of  indices  that  are  detected  as  nonredundant  at 
some  iteration  of  the  algorithm.  Thus  In  =  0  in  the  beginning  of  the  algorithm. 

We  denote  the  rational  matrix  in  (1.2 1  by  AT  and  the  scaled  matrix  defined  in 


(3.5)  by  AT .  The  matrix  B k  is  defined  by  (|4.1|)  and  is  composed  of  columns  a,j 
of  A,  where  j  £  I(xk,£),  while  the  matrix  Bn  is  composed  of  those  columns  of  A 
whose  indices  are  in  the  set  Thus  Bn  consists  of  the  normal  vectors  to  the 

nonredundant  constraints. 


Algorithm  for  constructing  the  set  of  generators  D k 

Let  the  current  iterate  Xk  and  a  parameter  e  >  0  be  given. 

%  Part  I.  Constructing  the  set  lN{xk,£ ) 


%  Constructing  the  working  index  set  I(xk,£ ) 
for  i  =  1  to  |  J| 

if  0  <bi  —  afxk  <  £ 

I(xk,£)  <-  v,  Bk+-  at; 

endif 

endfor 


if  rank  (Bk)  =  \I(xk,e)\ 

%  Bk  has  full  rank ,  hence,  all  constraints  are  nonredundant 
IN{xk,£)  <—  I{xk,e)]  Bn  <-  Bk ; 

else 

%  using  information  from  the  previous  iterations  of  the  algorithm 

IN{xk,e)  <-  {I(xk,£)f]IN}', 
if  {I(xk,e)\IN(xk,  e)}  ^  0 

%  there  are  active  constraints  that  have  not  been  identified 
%  at  the  previous  iterations 
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%  Identification  of  the  nonredundant  and  redundant  constraints 
for  j  =  1  to  \I(xk,£)\ 


%  the  first  strategy 

%  Pj(xk)  is  the  projection  of  Xk  onto  {a;  €  3?"  :  ajx  =  bj} 
Pj(xk)  =xk  +  o,j(bj  -  aj Xk)', 

if  dfPfixk)  <  hi  for  all  i  €  I\j 

%  equality  at  Pj{xk)  holds  for  only  the  jth  constraint 


IN{xk,e)  <-  j]  Bn  <—  ay, 

else  %  at  least  two  constraints  hold  with  equality  at  Pj{xk) 
%  or  the  point  Pj(xk )  is  infeasible ; 


%  the  second  strategy 


solve  LP  problem  (3.161;  let  x*  be  a  solution  to  (3.161; 


if  aj x*  <  bj  %  the  j  th  constraint  is  redundant 

take  ajX  <  bj  out  of  S!  and 
take  j  out  of  the  sets  I  and  I{xk,£)\ 
else  %  the  jth  constraint  is  nonredundant 
IN(xk,s)  <-  j ;  Bn  <-  aj] 

endif 


endif 

endfor 

endif 

endif 


%saving  information  for  the  next  iterations 
In  lN{xk,e)] 


%  Part  II.  Constructing  the  set  of  generators  Dk 
if  rank(BN)  =  |(IN)(xk,e)| 

%  nondegenerate  case 

[Dk]  =  Set.QR(BN)  or  [Dk]  =  SeCLU(BN ); 

else 

%  degenerate  nonredundant  case 
'  [Dk]  =  Set(BN)] 

endif 


Function  Set-QR  constructs  a  set  of  generators  by  using  the  QR  decomposition. 
The  procedure  is  defined  as  follows. 

function  [D]  =  Set_QR  (B) 

[Q,R]  =  qr(B ); 

r=rank(R); 

[Ql,  Q2,  f?l]  =  decomposition  (Q,  R,  r); 

D1  —  ±(Q1  *  R1  *  inv{RV  *  Rl))\ 

D2  —  ±Q2] 

D  =  [D1  D2}] 

end 
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Function  Set.LU  is  described  next  and  constructs  a  set  of  generators  by  using  the 
LU  factorization. 

function  [D]  =  Set  LU  (B) 

DL  *- ±(V  *  inv(V'  *  V)); 

D2  <-  ±(I  -  V  *  inv{V'  *V)  *  V'); 

D  =  [D1  D2]; 
end 


Function  Set  constructs  a  set  of  generators  in  the  case  when  the  matrix  B n  does 
not  have  full  rank.  If  Bn  has  tun  columns  and  rank  r,  then  the  number  S  of  the 
maximal  linearly  independent  subsets  of  the  columns  of  Bn  is  given  by  (4.5 1 . 

function  [D]  =  Set  (B) 
for  i  =  1  to  S 

Bi  is  composed  of  r  linearly  independent  columns  of  B ; 

[A]  =  Set-QR(Bi)  or  [A]  =  Set-LU(Bi)-, 


D  <-  D  n  A; 

%D  (~l  A  eliminates  the  identical  vectors  from  the  set  D  % 

endfor 

end 

To  reduce  the  number  of  operations  in  the  function  Set,  we  can  use  properties  of 
QR  factorizations.  Namely,  at  the  first  step,  a  QR  decomposition  of  B  is  computed. 
Then  at  the  ith  step,  i  =  2, ...  ,S,  only  one  of  the  first  r  columns  of  B  is  replaced 
with  some  j th  column  of  B,  where  j  >  r,  in  such  a  way  that  B ^  ^  Bi2  for  any  steps 
i\  and  i 2,  and  this  is  used  to  simplify  computing  another  QR  factorization  of  B.  The 
QR  factorization  of  B  is  used  to  construct  a  set  of  generators  for  A ,  since  A  can  be 
defined  as  the  first  r  columns  of  B  at  the  ith  step. 


In  the  next  two  tables,  we  present  some  numerical  results.  In  Table  4.1  we  report 
results  regarding  identification  of  the  nonredundant  and  redundant  constraints  at  the 
current  iterate  x k  and  regarding  constructing  the  set  of  the  nonredundant  constraints 
lN{xk,e).  In  Table  4.2  we  report  results  concerning  the  number  of  iterations  in 


the  GPS  algorithm  that  is  described  in  Fig|2.1|  and  is  applied  to  a  problem  with  an 
increasing  number  of  redundant  constraints. 


Table  4.1 

Constructing  the  set  I r^(x^,  e)  at  the  current  iterate 


variables 

\I{xk,e)\ 

e-active 

constraints 

|A(a;fc,e)| 

nonredundant 

constraints 

detected  as 
by  the  projection 
method 

nonredundant 
by  the  LP  program 

3 

6 

6 

6 

5 

6 

5 

4 

1 

5 

7 

7 

6 

1 

5 

7 

7 

5 

2 

In  Table  |4~2)  we  report  results  for  the  same  problem  with  no  redundant  constraints 
in  the  first  row,  with  one  additional  redundant  constraint  in  the  second  row,  and  with 
two  additional  redundant  constraints  in  the  third  row.  Since  the  algorithm  detects 
the  redundant  constraints  and  takes  them  out,  the  number  of  iterations  is  the  same 
in  all  three  tests. 
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Table  4.2 

GPS  algorithm  in  the  degenerate  case 


variables 

constraints 

redundant 

constraints 

iterations 

5 

9 

0 

18 

5 

10 

1 

18 

5 

11 

2 

18 

5.  Concluding  remarks.  It  is  interesting  to  compare  our  results  with  those  of 
Lewis  and  Torczon  [3]  and  of  Price  and  Coope  m-  Lewis  and  Torczon  [5j  proved  that 
it  is  possible  to  find  a  set  of  generators  for  the  cones  K°(xk,  e)  in  both  the  degenerate 
and  nondegenerate  case,  but  left  all  details  of  numerical  implementation  for  future 
work. 

Price  and  Coope  [T2]  presented  a  new  result  that  can  be  used  for  constructing  a 
set  of  generators  in  the  degenerate  case.  It  follows  from  their  result  that,  in  order  to 
construct  a  set  of  generators,  it  is  sufficient  to  consider  maximal  linearly  independent 
subsets  of  the  active  constraints.  However,  this  approach  implies  enumeration  of  all 
possible  linearly  independent  subsets  of  maximal  rank  and  does  not  take  into  account 
properties  of  the  problem  that  can  help  to  reduce  this  enumeration.  Price  and  Coope 
m  outlined  an  algorithm  for  constructing  frames,  but  did  not  consider  details  of  the 
numerical  implementation  in  the  degenerate  case. 

To  construct  the  set  of  generators,  we  first  classify  constraints  as  redundant  and 
nonredundant  by  applying  some  results  concerning  redundancy  from  the  mathemat¬ 
ical  programming  literature  and  by  using  our  approach  presented  in  the  paper.  We 
give  a  detailed  algorithm  for  constructing  the  set  of  generators.  However,  the  de¬ 
generate  nonredundant  case,  when  all  constraints  are  nonredundant  but  linearly  de¬ 
pendent,  still  implies  enumeration  of  all  maximal  linearly  independent  subsets  of  the 
constraints.  Therefore,  the  issue  left  for  future  analysis  is  whether  it  is  possible  to 
reduce  the  enumeration  of  the  subsets  in  the  degenerate  nonredundant  case. 

This  work  was  begun  at  the  IMA  while  the  first  author  was  a  postdoctoral  fellow 
and  the  second  was  a  longterm  visitor.  We  thank  the  IMA  for  providing  such  a  fine 
atmosphere  for  collaboration. 
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